Discontinuous Molecular Dynamics for Semi-Flexible and Rigid Bodies 
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A general framework for performing event-driven simulations of systems with semi-flexible or 
rigid bodies interacting under impulsive torques and forces is outlined. Two different approaches are 
presented. In the first, the dynamics and interaction rules are derived from Lagrangian mechanics in 
the presence of constraints. This approach is most suitable when the body is composed of relatively 
few point masses or is semi-flexible. In the second method, the equations of rigid bodies are used 
to derive explicit analytical expressions for the free evolution of arbitrary rigid molecules and to 
construct a simple scheme for computing interaction rules. Efficient algorithms for the search for 
the times of interaction events are designed in this context, and the handling of missed interaction 
events is discussed. 
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I. INTRODUCTION 

There has been an increasing interest over the last 
decade in performing large-scale simulations of colloidal 
systems, proteins, micelles and other biological assem- 
blies. Simulating such systems, and the phenomena that 
take place in them, typically requires a description of 
dynamical events that occur over a wide range of time 
scales. Nearly all simulations of such systems to date 
are based on following the microscopic time evolution 
of the system by integration of the classical equations 
of motion. Usually, due to the complexity of inter- 
molecular interactions, this integration is carried out in a 
step-by-step numerical fashion producing a time ordered 
set of phase-space points (a trajectory) . This information 
can then be used to calculate thermodynamic properties, 
structural functions or transport coefficients. An alter- 
native approach, which has been employed in many con- 
texts, is to use step potentials to approximate intermolec- 
ular interactions while affording the analytical solution of 
the dynamics [T]- [3- 0- EH- • The simplification in the in- 
teraction potential can lead to an increase in simulation 
efficiency since the demanding task of calculating forces 
is reduced to computing momentum exchanges between 
bodies at the instant of interaction. This approach is 
called event-driven or discontinuous molecular dynamics 
(DMD). 

In the DMD approach, various components of the sys- 
tem interact via discontinuous forces, leading to impul- 
sive forces that act at specific moments of time. As a 
result, the motion of particles is free of inter- molecular 
forces between impulsive events that alter the trajectory 
of bodies via discontinuous jumps in the momenta of the 
system at discrete interaction times. To determine the 
dynamics, the basic interaction rules of how the (linear 
and angular) momenta of the body are modified by col- 
lisions must be specified. 

For molecular systems with internal degrees of free- 



dom it is straightforward to design fully-flexible models 
with discontinuous potentials, but DMD simulations of 
such systems are often inefficient due to the relatively 
high frequency of internal motions 0. This inefficiency is 
reflected by the fact that most collision events executed 
in a DMD simulation correspond to intra rather than 
inter-molecular interactions. On the other hand, much 
of the physics relevant in large-scale simulations is insen- 
sitive to details of intra-molecular motion at long times. 
For this reason, methods of incorporating constraints into 
the dynamics of systems with continuous potentials have 
been developed that eliminate high frequency internal 
motion, and thus extend the time scales accessible to sim- 
ulation. Surprisingly, relatively little work has appeared 
in the literature on incorporating such constraints into 
DMD simulations. The goal of this paper is to extend 
the applicability of DMD methods to include constrained 
systems and to outline efficient methods that are gen- 
erally applicable in the simulations of semi-flexible and 
rigid bodies interacting via discontinuous potentials. 

In contrast to systems containing only simple spherical 
particles HIE IE 110, the application of DMD meth- 
ods to rigid-body systems is complicated by two main 
challenges. The first challenge is to analytically solve 
the dynamics of the system so that the position, veloc- 
ity, or angular velocity of any part of the system can 
be obtained exactly. This is in principle possible for a 
rigid body moving in the absence of forces and torques, 
even if it does not possess an axis of symmetry which 
facilitates its motion. However, an explicit solution suit- 
able for numerical implementation seems to be missing 
in the literature (although partial answers are abundant 
H EE IH IH El 0) For this reason, we will present 
the explicit solution here. Armed with a solution of the 
dynamics of all bodies in the system, one can calculate 
the collision times in an efficient manner, and in some 
instances, analytically. 

The second challenge is to determine how the impulsive 
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forces lead to discontinuous jumps in the momenta of the 
interacting bodies. For complicated rigid or semi-flexible 
bodies, the rules for computing the momentum jumps 
are not immediately obvious. It is clear however that 
these jumps in momenta must be consistent with basic 
conservation laws connected to symmetries of the under- 
lying Lagrangian characterizing the dynamics. Often the 
basic Lagrangian is invariant to time and space trans- 
lations, and rotations, and, hence, the rules governing 
collisions must explicitly obey energy, momentum, and 
angular momentum constraints. Such conservation laws 
can be utilized as a guide to derive the proper collision 
rules. 

A first attempt to introduce constraints into an 
event-driven system was carried out by Ciccotti and 
Kalibaeva|l5(, who studied a system of rigid, diatomic 
molecules (mimicking liquid nitrogen). Furthermore, 
non-spherical bodies of a special kind were treated by 
Donev et of. [161 Il7| by assuming that all rotational mo- 
tion in between interaction events was that of a spheri- 
cally symmetric body. More recently, a spherically sym- 
metric hard-sphere model with four tetrahedral short 
ranged (sticky) interactions (mimicking water) has been 
studied by De Michele et af.[18| with an event-driven 
molecular dynamics simulation method similar to the 
most basic scheme presented in this paper. This work 
primarily focuses on the phase diagram of this "sticky" 
water model as a prototype of network forming molecular 
systems. Our purpose, in contrast, is to discuss a gen- 
eral framework that allows one to carry out event-driven 
DMD simulations in the presence of constraints and, in 
particular, for fully general rigid bodies. The method- 
ology is applicable to modeling the correct dynamics of 
water molecules in aqueous solutions as well as other 
many body systems. 

The paper is organized as follows. Section ILTI discusses 
the equations of motions in the presence of constraints 
and Sec. II I II discusses the calculation and scheduling of 
collision times. The collision rules are derived in Sec. lIVI 
In Sec. E] it is shown how to sample the canonical and 
microcanonical ensembles and how to handle subtle is- 
sues concerning missing events that are particular to 
event-driven simulations. Finally, conclusions are pre- 
sented in Sec. IVII 



where the index a runs over all constraints present in the 
system and r n is a generalized vector whose components 
are the set of all Cartesian coordinates of the N total par- 
ticles in the system. For fully rigid bodies, the number 
of constraints c can easily be calculated by noting that 
the number of spatial degrees of freedom of an n-particle 
body is 3n in 3 dimensions, while only 6 degrees of free- 
dom are necessary to completely specify the position of 
all components of a rigid body: 3 degrees of freedom for 
the center of mass of the object and 3 degrees of freedom 
to specify its orientation with respect to some arbitrary 
fixed reference frame. There are therefore c = 3n — 6 
constraint equations for a single rigid body with n point 
masses. Below, these point masses will be referred to as 
atoms while the body as a whole will be called a molecule. 
A typical constraint equation fixes the distance between 
atoms i and j in the molecule to be some value, d, i.e., 
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The equations of motion for the system follow from 
Hamilton's principle of stationary action, which states 
that the motion of the system over a fixed time interval 
is such that the variation of the line integral defining the 
action S is zero: 



S 



C(r N (t),r N (t))dt 



(3) 



SS = 6 C(r N (t),r N (t))dt = 0, (4) 

where the Lagrangian C in the presence of the constraints 
is written in Cartesian coordinates as 

N 

C(r N ,r N ) = ^2 -^W 2 ~ ^ a o- a {r N ) - $(r N ), (5) 



where $ is the interaction potential. For clarity through- 
out this paper, the Einstein summation convention will 
be used for sums over repeated greek indices i.e. \ a a a = 
J2a=i ^aO~ a , whereas the sum over atom indices will be 
written explicitly. In Eq. ©, the parameters X a are La- 
grange multipliers to enforce the distance constraints o a . 
The resulting equations of motion are: 



II. EQUATIONS OF MOTION WITH 
CONSTRAINTS 



A. Constrained dynamics 

The motion of rigid bodies can be considered to be a 
special case of the dynamics of systems under a mini- 
mal set of c time-independent holonomic constraints (i.e. 
dependent only on positions) that fix all intra-body dis- 
tances: 



<7 a {r N ) = 0, 



(1) 



d dC 
dt dri 

ms.fi = 



dC 
dri 

da a d$ 
a dr t dri ' 



(6) 
(7) 



For an elementary discussion of constrained dynamics 
in the Lagrangian formulation of mechanics, we refer to 
Ref.lU 

When there are no interactions, such as for a single 
molecule, the potential <& = and Eq. Q) becomes 



m,r, 



-A, 



(8) 
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These equations of motion must be supplemented by 
equations for the c Lagrange multipliers X a , which are 
functions of time. Although the X a are not functions of 
t*jv and r n in a mathematical sense, it will be shown 
below that once the equations are solved they can be ex- 
pressed in terms of r^v and rjv. Note that the equations 
of motion show that even in the absence of an external 
potential, the motion of the point masses (atoms) mak- 
ing up a rigid body (molecule) are non-trivial due to the 
emergence of a constraint force —X a da a / 'dvi. 

In fortuitous cases, the time dependence of the La- 
grange multipliers is relatively simple and can be solved 
for by Taylor expansion of the Lagrange multipliers in 
time t. To evaluate the time derivatives of the multipli- 
ers, one can use time derivatives of the initial constraint 
conditions, which must vanish to all orders. The result 
is a hierarchy of equations, which, at order k, is linear in 

(k) 

the unknown fcth time derivatives Aq but depends on 
the lower order time derivatives Xa^ , Xa , ■ ■ ■ X^ ^ . In 
exceptional circumstances, this hierarchy naturally trun- 
cates. For example, for a rigid diatomic molecule with a 
single bond length constraint, one finds that the hierar- 
chy truncates at order k = 0, and the Lagrange multiplier 
is a constant However this is not the typical case. 

Alternatively, since the constraints o~ a (y*jv) — are to 
be satisfied at all times t, and not just at time zero, their 
time derivatives are zero at all times. From the first time 
derivative 

o a {r N ) = 0, 

one sees that the initial velocities Vi — must obey 



Edo a 



(9) 



for each constraint condition a. The Lagrange multipli- 
ers can be determined by the condition that the second 
derivatives of all the constraints vanish so that 



a- a (r N ) = 



da a (r N ) 



dr. 



d 2 a a (r N ) . 

rt + > ta ■ — — n 



drjdri 



E 



da a (r N ) A/? dap(r N ) 



dn mi dr l 
d 2 a a {r N ) . 



drjdri 



= 0, 



yielding a linear equation for the Lagrange multipliers 
that can be solved in matrix form as 



X a (t) = Z-*(r N (t))Tp(r N (t),r N (t)), (10) 



where 



T (r Nl r N ) = 



d 2 a l3 {r N ) 
drjdri 



Z afs {r N ) = ^ 



1 da a (r N ) do-p(r N ) 



m, dr,. 



dr t 



(11) 



(12) 



It may be shown that with A Q given by Eq. IjlOfl. all 
higher time derivatives of o~ a are automatically zero. 

As Eq. (|10(1 shows, in general the Lagrange multipliers 
are dependent on both the positions r^ and the velocities 
rpf of the particles. To see that this makes the dynamics 
non-Hamiltonian, the equations of motion can be cast 
into Hamiltonian-like form using pi — rriiri, i.e., 

Pi 



m, 



Pi 



-X a 



dag 

dn ' 



(13) 



where it is apparent that the forces in the system de- 
pend on the momentum through A in Eq. (|1C)|> . There 
exists no Hamiltonian that generates these equations of 
motion. |2l| 

Since the underlying dynamics of the system is 
non-Hamiltonian, the statistical mechanics of the con- 
strained system is potentially more complex. In general, 
phase-space averages have to be defined with respect to 
a metric that is invariant to the standard measure of 
Hamiltonian systems, but dr^dpN is not conserved un- 
der the dynamics and the standard form of the Liouville 
equation does not hold [22L Vl?\ . In general, there is a 
phase-space compressibility factor k associated with the 
lack of conservation of the measure that is given by mi- 
nus the divergence of the flow in phase space. It may be 
shown that pa. 12^ 

*=-^-|r^ = l ln||z(rjv)l1 ' 

where ||Z(rjv)|| is the determinant of the matrix Z Q ^(r7v) 
defined in Eq. H12|) . The compressibility factor is 
related to the invariant pha se-space metric d\x = 
y / g(r N ,p N ,t)dr N dp N withjH HI 



Vg( r N,PN,t) = ||Z(rjv)||. 



(14) 



Statistical averages are therefore defined for the 
non-Hamiltonian system as |2rjl27| 

(X(r N ,p N )) = i J dp N dr N ||Z(rjv)|| 
x X(r N ,p N ) p(r N ,p N ) 
x Y\ {r N ))S(a a (r N , p N )) ,(15) 

ot 

where p(r^,p^) is the probability density for the uncon- 
strained system and Q is the partition function for the 
constrained system, given by 



Q 



dp N dr N ||Z(rjv)|| p(r N ,PN) 
Y\_S(a a (r N ))S {<J a (r N ,p N )) . 



Although the invariant metric is non-uniform for many 
constrained systems, for entirely rigid systems the Z(r^) 
matrix is a function only of the point masses and fixed 
distances. Hence the term ||Z(rjv)|| acts as a multiplica- 
tive factor which cancels in the averaging process. 
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B. Free motion of rigid bodies 

Although the solution of the dynamics of constrained 
systems via time-independent holonomic constraints is 
intellectually appealing and useful in developing a formal 
statistical mechanics for these systems, it is often difficult 
to analytically solve for the values of the Lagrange mul- 
tipliers at arbitrary times. One therefore often resorts to 
numerical solutions of the multipliers in iterative form, 
using algorithms such as SHAKE 28]. Such an approach 
is not really consistent with the principles of DMD, in 
which a computationally efficient means of calculating 
event times is one of the great advantages of the method. 
For fully-constrained, rigid bodies, it is more sensible to 
apply other, equivalent, approaches, such as the principal 
axis or quaternion methods, to calculate analytically the 
evolution of the system in the absence of external forces. 

The basic simplification in the dynamics of rigid bodies 
results from the fact that the general motion of a rigid 
body can be decomposed into a translation of the center 
of mass of the body plus a rotation about the center of 
mass. The orientation of the body relative to its center of 
mass is described by the relation between the so-called 
body frame, in which a set of axes are held fixed with 
the body as it moves, and the fixed external laboratory 
frame. The two frames of reference can be connected by 
an orthogonal transformation, such that the position of 
an atom i in a rigid body can be written at an arbitrary 
time t as: 

r i (t) = R(t)+A^(t)f i) (16) 

where r, is the position of atom i in the body frame 
(which is independent of time), R(t) is the center of mass, 



and the matrix At (t) is the orthogonal matrix that con- 
verts coordinates in the body frame to the laboratory 
frame. Note that matrix-vector and matrix-matrix mul- 
tiplication will be implied throughout the paper. The 
matrix A^(t) is the transpose of A(t), which converts co- 
ordinates from the laboratory frame to the body frame at 
time t. The elements composing the columns of the ma- 
trix A^ (t) are simply the coordinates of the axes in the 
body frame written in the laboratory frame. Note that 
Eq. I jl 6|l implies that the relative vector r , (t) satisfies 

r i =r i -R = A*f i , (17) 

Here as well as below, we have dropped the explicit time 
dependence for most time dependent quantities with the 
exception of quantities at time zero or at a time that is 
integrated over. 

One sees that in order to determine the location of 
different parts of the body in the laboratory frame, the 
rotation matrix A must be specified. This matrix satis- 
fies a differential equation that will now be derived and 
subsequently solved. 

Before doing so, it will be useful to restate some prop- 
erties of rotation matrices and establish some notation 
to be used below. Formally, a rotation matrix U is an 
orthogonal matrix with determinant one and whose its 
inverse is equal to its transpose TJ'. Any rotation can 
be specified by a rotation axis n = (711,712,713) and an 
angle 9 over which to rotate. Here h is a unit vector, so 
that one may also say that any non-unit vector 9n can 
be used to specify a rotation, where its norm is equal to 
the angle and its direction is equal to the axis n. Ac- 
cording to Rodrigues' formula, the matrix corresponding 
to this rotation is[20| 



n \ + ( n 2 + n i) cos ^ Hi7i2 (1 — cos 9) — 77.3 sin 9 711713 (1 — cos 9) + 112 sin 9^ 
TJ (On) — \ n\ri2 (1 — cos 0) + n 3 sin 9 n 2 + (nf + n|) cos n 2 n 3 (1 — cos 9) — rii sin ] . (18) 
yn^rii (1 — cos 0) — n 2 sin 9 n 3 n 2 (1 — cos 9) + n\ sin 9 n\ + (n\ + n|) cos 9 



The derivation of the differential equation for A starts 
by taking the time derivative of Eq. (JTSJ, yielding 

Vi - V = AJ fi = A f Ar,. (19) 

From elementary classical mechanics |20j. it is known that 
this relative velocity can also be written as 

Vi-V = 0jxri, (20) 

where uj is the angular velocity vector in the lab frame. 
Since both Eq. I|19fl and Eq. I|2l)|l are true for any vector 
r, , it follows that At A is the matrix representation of a 



cross product with the angular velocity lj, i.e., 

(0 -u> 3 u 2 \ 
lo 3 -W! = W(w). (21) 
LUl J 

Multiplying Eq. I|21|) on the right with At and taking the 
transpose on both sides (note that W is antisymmetric) 
yields 

A = -AW(w). (22) 

This equation involves the angular velocity in the labo- 
ratory frame, but the rotational equations of motion are 
more easily solved in the body frame. The angular ve- 
locity vector transforms to the body frame according to 
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Cj = A< 



(23) 



For any rotation A and vector x one has W(Ax) 
AW(cc) A T , hence one can write 



W(u>) = W(A f w) = A f W(w) A. 



(24) 



Substituting Eq. 
equation for A: 



into Eq. (|22|l yields the differential 



-W(u>) A. 



(25) 



Although the choice of body frame is arbitrary, per- 
haps the most convenient choice of axes for the body are 
the so-called principal axes in which the moment of iner- 
tia tensor I is diagonal, i.e., I = diag(/i, I 2 , I 3 ). Choosing 
this reference frame as the body frame, the representa- 
tion of the components of the angular momentum L is 



(26) 



where Ik and u>k are the principal moments of inertia and 
principal components of the angular velocity. 

The time dependence of the principal components of 
the angular velocity may be obtained from the standard 
expression for the torque in the laboratory frame: 




t = L = A+L + Ati 

= A f u x L + L 



(27) 



where Eq. I|25|l was used in the last equality. Trans- 
forming Eq. (|27|l to the principal axis frame gives Euler's 
equations of motion for a rigid body 



I 2 U) 2 - ^1^3(^3 - h) = T 2 

I3W3 - u!iu> 2 (h - h) = f 3 , 



(28) 



where Tfe are the components of the torque f = At in 
the body frame. Note that even in the absence of any 
torque, the principal components of the angular velocity 
are in general time dependent. 

Once the angular velocity u) is known, it can be sub- 
stituted into Eq. (|25|l for the matrix A. The general 
solution of Eq. (|25|l is of the form 



A = PA(0). 



(29) 



where P is a rotation matrix itself which 'propagates' the 
orientation A(0) to the orientation at time t. P satisfies 
the same equation (|25|) as A, but with initial condition 
P(0) = 1. By integrating this equation, one can obtain 
an expression for P. At first glance, it may seem that 
P can only be written as a formal expression containing 



a time-ordered exponential. However, for the torque-free 
case r = 0, the conservation of angular momentum and 
energy and the orthogonality of the matrix P can be used 
to derive the following explicit expression|29j (implicitly 
also found in Ref. WS): 



TiT 2 



(30) 



Here Ti and T 2 are two rotation matrices. The matrix 
Ti rotates L(Q) to L and can be written as 



where 



cic 2 - sis 2 c 3 

-SlC 2 - C!S 2 C 3 

S2S3 



C1S2 + S1C2C3 

SiS 2 + C X C 2 C 3 

-c 2 s 3 



S1S3 
C1S3 
C3 



si 



s 2 



S3 



h. 

Li(0) 
L X (0) 
L ± L 3 (0) 



L 3 L ± (0) 



V- 



c 2 = 



C3 = 



L± 
L 2 (0) 



(31) 



(32) 



(33) 



L ± (0) 
L ± L ± (0)+L 3 L 3 (0) 



L' 2 



(34) 



and L± = y L\ + L\ while L = \L\. 

The matrix T2 can be expressed, using the notation in 
Eq. HH, as 



U(-^ 1 i(0)), 



where the angle ip is given by 



i<= / dt'tt{t'), 
Jo 



with 



n = l 



hu>i + I 2 u> 2 



(35) 



(36) 



(37) 



The angle tp can be interpreted as an angle over which the 
body rotates. If the body rotates one way, the laboratory 
frame as seen from the body frame rotates in the opposite 
way, which explains the minus sign in Eq. l| 35l) . For the 
derivation of Eqs. I|3U|I - (|37|I we refer to Ref. |2JJ- Similar 
equations, but in a special reference frame, can be found 
in Ref. Q 

In the following, the solution of Eq. 128fl with r = 
for bodies of differing degrees of symmetry will be ana- 
lyzed and then used to obtain explicit expressions for the 
matrix P as a function of time and of the initial angular 
velocity in the body frame w(Q). 



Case 1. Spherical rotor 

For the case of a spherical rotor in which all three 
moments of inertia are equal, I\ = I 2 = I3, the form of 
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the Euler equations l|28[l is particularly simple: IiQj = 0. 
It is therefore clear that all components of the angular 
velocity in the body frame are conserved, as are those 
of the angular momentum. As a result, Ti in Eq. (|3 1|) 
is equal to the identity matrix. A second consequence is 
that Jl in Eq. I|37|l is constant, so that ip = fit where il 
may be rewritten, using I\ = I2 — I3, as SI = L/1% = \u)\. 
Therefore Eqs. and give 



U(-wt) 



(38) 



corresponding to a rotation by an angle of —Clt around 
the axis w/O. 



Case 2. Symmetric top 

For the case of a symmetric top for which I\ = I2, 
one can solve the Euler equations (|28[) in terms of simple 
sines and cosines, since Eq. (|28|l becomes 



Cj 2 = -UJ p LOi 
LU 3 = 0, 



(39) 



where oj p = ^1 — -^-jil^O) is the precession frequency. 
The full solution of the Euler equations 139fl is given by 

u>i (0) cos ujpt + a>2 (0) sin ui p t 
u> = l — o>i(0) sin Upt + 0)2(0) cos tu p t | . (40) 
£ 3 (0) 

Using Eq. (|40|) and the fact that L± and L3 are conserved 
in this case, one can easily show that Ti is given by 

cos LOpt sin ujpt 0\ 
Ti = I — sincjpt cosuipt = U(— u v tz). (41) 
' l) 

and one can determine $1 from Eq. (|37|l : 

_ L[hul(0) + / 1 £j(0)] _ L 

This is a constant so that il) = ^-t. Thus 

L(0)t\ 



T 2 = U 



and one gets from Eq. (|30J) : 



h 



(43) 



P = U(-^)u(-^ 



(44) 



Case 3. Asymmetric body 

If all the principal moments of inertia are distinct, the 
time dependence of the angular velocity u) involves ellip- 
tic functions [23]. While this may seem complicated, effi- 
cient standard numerical routines exist to evaluate these 
functions [30L l3lL I32I l33l l3^ | . More challenging is the eval- 
uation of the matrix P. While its exact solution has been 
known for more than 170 yea rs El ITfij. it is formulated— 
even in more recent texts [ill Il2j — in terms of undeter- 
mined constants and using complex algebra, which hin- 
ders its straightforward implementation in a numerical 
simulation. It is surprisingly difficult to find an explicit 
formula in the literature for the matrix P as a function of 
the initial conditions, which is the form needed in DMD 
simulations. For this reason, the explicit general solution 
for P will briefly be presented here in terms of general 
initial conditions. The details of the derivation can be 
found elsewhere |29j. 

Following Jacobi^3> ^ is useful to adopt the conven- 
tion that I2 is the moment of inertia intermediate in mag- 
nitude (i.e., either I\ < 1% < I3 or I\ > I 2 > I3) and one 
chooses the overall ordering of magnitudes, such that: 



h > h > h if E R > — 
Z12 



h<h< h if E R < 



2h 



(45) 



where En is the rotational kinetic energy Er = ^{IiCj{ + 
ui\ + I3L1J3) and L is the norm of the angular momentum 
L = (IiLof + + ifd)!) 1 / 2 . Without this convention 
some quantities defined below would be complex valued, 
which is numerically inconvenient and inefficient. Note 
that in a simulation molecules will often be assigned a 
specific set of physical inertial moments with fixed order, 
i.e. not depending on the particular values of Er and 
L. A simple way to nevertheless adopt the convention in 
Eq. 145|) is to introduce internal variables a> int = Tj Tlt o>, 
I mt = T int IT int and A int = T mt AT mt which differ 
when necessary from the physical ones by a rotation given 
by the rotation matrix 



(46) 



This matrix interchanges the x and z directions and re- 
versed the y direction, and is equal to its inverse. 

The Euler equations l|28|l can be solved because there 
are two conserved quantities E and L 2 which allow Cj\ 
and U3 to be expressed in terms of il>2, at least up to 
a sign which the quadratic conserved quantities cannot 
prescribe. In this way the three coupled equations lt2~8|) 
are reduced to a single ordinary differential equation for 
du)2/dt 1 from which t can be solved as an integral over 
H12 ■ This is an incomplete elliptic integral of the first 
kind |30l l3l| . To get u>2 as a function of t, one needs its 
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0.5- 



-0.5 




cn(x;k) 

- — sh(a;A') 

■ ■ ■ dn(x;k) 

cos(7l^/(2JO) 

sin(7l;t/(2X)) 



FIG. 1: Example of the elliptic functions cn (solid line), sn 
(bold dashed line), dn (dotted line) for m = 0.465 (K = 2.05, 
K' — 1.72, q — 0.071). Also plotted are the cosine (short 
dashed line) and sine (thin short dashed line) with the same 
period, for comparison. 



inverse, which is the elliptic function sn|3Jj,|31(. Without 
giving further details, the solution of the Euler equations 
is given by[nim[ll,|2| 



(Jim cn(oj p t + e\m) 
(j 2m sii(uj p t + e\m) 
, <j 3m dn(ui p t + e\m) 



(47) 



Here cn and dn are also elliptic functions [3(J, Kill l35| . 
while the Ui m are the extreme (maximum or minimum) 
values of the u>i and are given by 



(Jim = Sgn(d>i(0))4 



U2r, 



sgn(J)i(0)) 1 



u 3m = sgn(w 3 (0)) 1 



< L 2 -2I 3 E R 
h(h-h) 

/L 2 - 2I 3 E R 
Hh-h) 

! L 2 - 2hE R 
h(h-h) ' 



(48) 



where sgn(x) is the sign of x. Furthermore, in Eq. 147|) 
the precession frequency lo v is given by 



uj p = sgn(/ 2 - ^3) sgn(a) 3 (0))^ 



(L 2 - 2hE R ){I 3 - h 



hhh 



The elliptic functions are periodic functions of their first 
argument, and look very similar to the sine, cosine and 
constant function. They furthermore depend on the el- 
liptic parameter m (or elliptic modulus y/m), which de- 
termines how closely the elliptic functions resemble their 
trigonometric counterparts, and which is given by 



(h - I 2 )(L 2 - 2I 3 E R ) 
{I 3 -h){L^-2IiE R y 



(50) 



By matching the values of 012 at time zero, one can 
determine the integration constant e: 



F(uj 2 o/oj2m\m), 



(51) 



where F is the incomplete elliptic integral of the first 
kindlMISl 



F{y\m) 



dx 



y/ (1 — x 2 )(l — mi 2 ) 



(52) 



In fact, sn(x|m) is simply the inverse of this function. 
As a result of the ordering convention in Eq. I|45|l . the 
parameter m in Eq. (|50|l is guaranteed to be less than 
one, which is required in order that F{y\m) in Eq. I|52|) 
not be complex-valued. 

Three more numbers can be derived from the elliptic 
parameter m which play an important role in the elliptic 
functions. These are the quarter-period K — F(l\m), 
the complementary quarter-period K' — — m ) an d 

the nome q — exp(—TrK'/K), which is the parameter in 
various series expansions. 

The period of the elliptic functions cn and sn is equal 
to 4K, while that of dn is 2K. These elliptic functions 
have the following Fourier series [30ll3l|: 



cn(x\m) 
sn(x\m) 
dn(cc|m) 



2tt ^ q n+1 ' 2 (2n+l)nx 



Tl = 



2K 



2ir ^ q n+1 ' 2 . (2n + l)7ra; 

?^ sm — (54) 



n=0 
7T 27T 

2K + ~K 



2K 



St 



q nnx 
s— cos . 

-q 2n K 



(55) 



Note that the right-hand side of Eqs. (|53|l - 155() depends 
on m through K = F(\\m) and q = exp[— 7rF(l|l — 
m)/K]. For m = 0, one gets q = and cn, sn and dn 
reduce to cos, sin and 1, respectively. The constancy of 
dn(x|m = 0) is reminiscent of the conservation of Cj 3 in 
the case of the symmetric top, and, indeed, for I\ = h, 
m = according to Eq. I|5(J|) . 

Typical values for q are quite small, hence often the 
elliptic function sn, cn and dn resemble the cos, sin and 
a constant function with value one (as e.g. in Fig.^J. F° r 
small values of q, the series expressions for the elliptic 
functions converge quickly (although this is not the best 
way to compute the elliptic functions j30l l3ll 1321 133 . l34| ) . 

Having given the solutions of the Euler equations, 
we now turn to the solution of Eq. (|25(l as given by 
Eqs. H29fl - H37fl . The expression on the right-hand side of 
Eq. (|37|l isnot a constant in this case but involves ellip- 
tic functions. Despite this difficulty, the integral can still 
be performed using some properties of elliptic functions, 
with the result 291 



i/j(t) =Ai+A 2 t-4>(t). 



(56) 



The constants Ai , A2 and the periodic functi on 4>(t) can 
be expressed using the theta function H(u\m) [Tot l3Cl l31| 
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<f>(t) = arg H(u)pt + e — irj\m) 
A\ = (f)(0) = argiJ(e — ir]\m) 
L dlogH(in\m) 



A, = — 



where we have used the definition 



i] = sgn(a> 30 )if' - F 



L 



(57) 
(58) 

(59) 



(60) 



The equations l|57|) - (|59(l involve complex values which 
are not convenient for numerical evaluation. Using the 
known series expansions of the theta function H and its 
logarithmic derivative [30j, l31| in terms of the nome q, 
these equations may be rewritten in a purely real form. 
In fact, one readily obtains the sine and cosine of "0, which 
are all that is needed in Eqs. I|18|) and iJSSJ, 



cos ip(t) 
sin ip(t) 



h r (t) cos(Ai + A 2 t) + hi(t) sm(A x + A 2 t) 

7WW+W) 

h r (t) sin(Ai + A 2 t) - hi(t) cos(A 1 + A 2 t) 



(61) 



with 



h r (t) — Re H(u)pt + e — \rj\m) 

oc 

= 2 g 1 / 4 ^(-l) l y i (" +1) cosh 



(2n + 1)7117 



n=0 



x sm 



2K 

(2n + l)n(w p t + e) 



2K 



(62) 



hi(t) = 



(2n + l)irr) 



Im H (u)pt + e — iry|m) 

oo 

-2 (Z 1 /4^ ( _ 1)> n(n + D sinh ^ 

(2n + l)w(u p t + e) 



n=0 



X COS ■ 



2K 



while the constant A\ is 

A\ = arctan[/ii(0)//i r (0)] + mr, 



(63) 



(64) 



where n = if h r (0) > 0, n = 1 if h r (0) < and /ij(0) > 
0, and n = -1 if h r (0) < and hi(0) < 0. Finally the 
constant A 2 is given bv [30ll3l| 



Ao 



L 
h 



IK UJp 

~2K 



j+1 



i 



,2n 



(65) 



where ^ = exp^^/if ). The series expansion in q in 
Eq. (|65|l convergences for £q 2 < 1. Because —if < rj < 
K' (cf. Eq. ||BUJ|), one has ^g 2 < g < 1, and the series 
always converges. Since q is typically small, the con- 
vergence is rarely very slow (e.g. for convergence up to 
relative order 8 one needs O (log 8/ log q) terms). Note 



that since the constants A\ and A 2 depend only on the 
initial angular velocities, they only need to be calculated 
once at the beginning of the motion of a free rigid body. 
On the other hand, the series expansions in Eqs. I|62|l and 
(|63|l . which have to be evaluated any time the positions 
are desired, have extremely fast convergence due to the 
gn(n+i) a pp ear j n g m these expressions (for example, un- 
less m > 0.95, the series converges up to O(10~ 15 ) occurs 
taking only three terms). 

There are efficient routines to calculate the functions 
cn, sn dn and F, see e.g. Refs. I30L l3ll l32l l33t l34l and the 
series in Eqs. (|o^|l and (|u3|) converge, the former 

two quite rapidly in fact. Therefore, despite an apparent 
preference in the literature for conventional numerical in- 
tegration of the equations of motion via many successive 
small time steps even for torque-free cases, the analytical 
solution can be used to calculate the same quantities in 
a computationally more efficient manner requiring only 
the evaluation of special functions. The gain in efficiency 
should be especially pronounced in applications in which 
many evaluations at various times could be needed, such 
as in the root searches in discontinuous molecular dy- 
namics (see below). 



III. DETECTION OF INTERACTION EVENTS 

A. Calculation of interaction events 

If the interaction potential between atoms i and j is 
assumed to be discontinuous, say of the form 



U if|n-r,|<d 
JU $! if \r i -r j \>d, 



(66) 



then rigid molecules interacting via this potential evolve 
freely until there is a change in the potential energy of 
the system and an interaction event or collision event 
occurs. The time at which an event occurs is governed 
by a collision indicator function = \rj — ri\ 2 — d 2 
defined such that at time t c , fij(t c ) = 0. Here, the time 
dependence of rj and can be obtained using the results 
of Sec. HD 

The simplest example of this kind of system consists 
of two hard spheres of diameter d located at positions rj 
and r.;. If two spheres are approaching, when they get to 
a distance d from one another, the potential energy would 
change from $i to $o = oo if they kept approaching one 
another. As this eventually would lead to a violation of 
energy conservation, the spheres bounce off one another 
in a hard-core collision at time t c , where t c is determined 
by the criterion f tj (t c ) = \rj(t c ) - n(i c )| 2 - d 2 = 0, i.e. 
by the zeros of the collision indicator function. Another 
kind of interaction event, with $o and $i finite, will be 
called a square-well collision because the potential then 
has a square well shape. 

To find the times at which collisions take place, the 
zeros of the collision indicator functions must be deter- 
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mined, which generally has to be done numerically. The 
calculation of the collision times of non-penetrating rigid 
objects is an important aspect of manipulating robotic 
bodies, and is also an important element of creating real- 
istic animation sequences. As a result, many algorithms 
have been proposed in these contexts to facilitate the 
event time search [36ll3^ . 

The search for the earliest collision event time can 
be facilitated using screenin g st rategies to decide when 
rigid bodies may overlap pa 137]. Usually, these involve 
placing the bodies in bounding boxes and using an effi- 
cient method to determine when bounding boxes inter- 
sect. The simplest way to do this in a simulation of rigid 
molecules is to place each molecule in the smallest sphere 
around its center of mass containing all components of 
the molecule |l5|. The position of the sphere is deter- 
mined by the motion of the center of mass, while any 
change in orientation of the rigid molecule occurs within 
the sphere. Collisions between rigid molecules can there- 
fore only occur when their encompassing spheres overlap, 
and the time at which this occurs can be calculated an- 
alytically for any pair of molecules. This time serves as 
a useful point to begin a more detailed search for colli- 
sion events (see below). Similarly, one can also calculate 
the time at which the spheres no longer overlap, and use 
these event times to bracket a possible root of the col- 
lision indicator function. It is crucial to make the time 
bracketing as tight as possible in any implementation of 
DMD with numerical root searches because the length 
of the time bracketing interval determines the required 
number of evaluations of the positions and velocities of 
the atoms, and therefore plays a significant role in the 
efficiency of the overall procedure. 

The simplest reliable and reasonably efficient means of 
detecting a root is to perform a grid search that looks for 
changes in sign of fij, i.e., one looks at fij(t + nAt) and 
fij (t+(n+l)At) for successive n. The time points t+nAt 
will be called the grid points. When a time interval in 
between two grid points is found in which a sign change 
of fij occurs, the Newton-Raphson algorithm 32J can be 
called to numerically determine the root with arbitrary 
accuracy. Since the Newton-Raphson method requires 
the calculation of first time derivatives, one must also 
calculate, for anytime t, the derivative dfij/dt = 2rij-Vij, 
where the notation m = rj — r, and «y = Vj — Vi has 
been used. Such time derivatives are readily evaluated 
using Eqs. JTJJ and l(T5)l, 

Unfortunately, while the Newton-Raphson method is a 
very efficient algorithm for finding roots, it can be some- 
what unstable when one is searching for the roots of an 
oscillatory function. For translating and rotating rigid 
molecules, the collision indicator function is indeed oscil- 
latory due to the periodic motion of the relative orienta- 
tion of two colliding bodies. It is particularly easy to miss 
so-called grazing collisions when the grid search interval 
At is too large, in which case the indicator function is 
positive in two consecutive points of the grid search, yet 
nonetheless "dips" below zero in the grid interval. It is 



important that no roots are missed, for a missed root can 
lead to a different, even infinite energy (but see Sec. |V] 
below) . To reduce the frequency of missing grazing colli- 
sions to zero, a vanishingly small grid interval At would 
be required. Of course such a scheme is not practical, 
and one must balance the likelihood of missing events 
with practical considerations since several collision indi- 
cator functions need to be evaluated at each point of the 
grid. Clearly the efficiency of the root search algorithm 
significantly depends on the magnitude of grid interval. 

To save computation time, a coarser grid can be uti- 
lized if a means of handling grazing collisions is im- 
plemented. Since the collision indicator function has a 
local extremum (maximum or minimum, depending on 
whether |r.; — rj\ 2 is initially smaller or larger than d 2 ) 
at some time near the time of a grazing collision, a rea- 
sonable strategy to find these kind of collision events is to 
determine the extremum of the indicator function in cases 
in which the indicator function fij itself does not change 
sign on the interval but its derivative dfij/dt does. Fur- 
thermore, since the indicator function at the grid points 
near a grazing collision is typically small, it is fruitful to 
search only for extrema when the indicator function at 
one of the grid points lies below some threshold value [38|. 
To find the local extrema of the indicator function, any 
simple routine of locating the extrema of a non-linear 
function can be utilized. For example, Brent's minimiza- 
tion method 32, 39], which is based on a parabolic inter- 
polation of the function, is a good choice for sufficiently 
smooth one-dimensional functions. Once the extremum 
is found, it is a simple matter to decide whether or not a 
real collision exists by checking the sign of fij. 

Once the root has been bracketed (either through 
a sign change of fij during the grid search or after 
searching for an extremum), one can simply use the 
Newton-Raphson algorithm to find the root to desired 
accuracy, typically within only a few iterations. The time 
value returned by the Newton-Raphson routine needs to 
be in the bracketed interval and dfij j dt < if /y was ini- 
tially positive and dfij/dt > if it was initially negative. 
If those criteria are not satisfied, the Newton-Raphson 
algorithm has clearly failed and a less efficient but more 
reliable method is needed to track down the root. For ex- 
ample, the Van Wijngaarden-Dekker-Brent method|32l 
[3£j > which combines bisection and quadratic interpola- 
tion, is guaranteed to converge if the function is known 
to have a root in the interval under analysis. 



B. Scheduling events 

In the previous section is was shown how to determine 
the time Uj at which two atoms i and j collide under the 
assumption that there is no other earlier collision. This 
we will call a possible collision event. In a DMD simula- 
tion, once the possible collision events at times have 
been computed for all possible collision pairs i and j, the 
earliest event t*, ■» = min^- should be selected. After 
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the collision event between atoms i* and j* has been exe- 
cuted (according to the rules derived in the next section), 
the next earliest collision should be performed. However, 
because the velocities of atoms of the molecules involved 
in the collision have changed, the previously computed 
collision times involving these molecules arc no longer 
valid. The next event in the sequence can be determined 
and performed only after these collision times have been 
recomputed. 

This process describes the basic strategy of DMD, 
which without further improvements would be needlessly 
inefficient. For if M is the number of possible collision 
events, finding the earliest time would require O(M) 
checks, and M — 0(N 2 ), while the number of invali- 
dated collisions that have to be recomputed after each 
collision would be O(N). Since the number of collisions 
in the system per unit of physical time also grows with 
N, the cost of a simulation for a given physical time 
would be 0(N 2 ) for the computation of collision times 
and 0(N 3 ) for finding the first collision event [icj. Fortu- 
nately, there are ways to significantly reduce this compu- 
tational cost 0, Ell E2, 113, E3, |4^| . The first technique, 
also used in molecular dynamics simulations of systems 
interacting with continuous potentials, reduces the num- 
ber of possible collision times that have to be computed 
by employing a cell division of the systempj. Note that 
while the times of certain interaction events (e.g. involv- 
ing only the molecule's center of mass) can be expressed 
in analytical form and thus computed very efficiently, 
the atom-atom interactions have, in general, an orien- 
tational dependence and the possible collision time has 
to be found by means of a numerical root search as ex- 
plained in the previous section. As a consequence, the 
most time consuming task in a DMD simulation with 
rigid bodies is the numerical root search for the collision 
times. One can however minimize the required number of 
collision time computations by dividing the system into 
a cell structure and sorting all molecules into these cells 
according to the positions of their centers of mass. Each 
cell has a diameter of at least the largest "interaction 
diameter" of a molecule as measured from its center of 
mass. As a result, molecules can only collide if they are 
in the same cell or in an adjacent cell, so the number 
of collision events to determine and to recompute after a 
collision is much smaller. In this technique, the sorting of 
molecules into cells is only done initially, while the sort- 
ing is dynamically updated by introducing a cell- crossing 
event for each molecule that is also stored0,E3. Since 
the center of mass of a molecule performs linear motion 
between collision events, one can express its cell-crossing 
time analytically and therefore the numerical computa- 
tion of that time is very fast. 

The second technique reduces the cost of finding the 
earliest event time. It consists of storing possible colli- 
sion and cell-crossing events in a time-ordered structure 
called a binary tree. For details we refer to Refsja and 
l42l (alternative event scheduling algorithms exist [jj, liif 
but it is not clear which technique is generally the most 



efficient H3.) 

Finally, a third standard technique is to update the 
molecules' positions and velocities only at collisions (and 
possibly upon their crossing the periodic boundaries), 
while storing the time of their last collision as a prop- 
erty of the molecule called its local clock^^. Whenever 
needed, the positions and velocities at later times can 
be determined from the exact solution of force-free and 
torque-free motion of the previous Sec. [n] 

The use of cell divisions, a binary event tree to man- 
age the events, and local clocks is a standard practice in 
DMD simulations and largely improves the simulation's 
efficiency . To see this, note that in each step of the sim- 
ulation one picks the earliest event from the tree, which 
scales as 0(log N) for randomly balanced trees [HE2- If 
it is a cell-collision event, it is then performed and subse- 
quently O(l) collisions and cell crossings are recomputed 
and added to the event tree (oc O (log N)). If it is a cross- 
ing event, the corresponding molecule is put in its new 
cell, new possible collision and crossing events are com- 
puted ((0(1)) and added to the tree (O (log N)). Then 
the program progresses to the next event. Since still 
O(N) real events take place per unit of physical time, 
one sees that using these techniques, the computational 
cost per unit of physical time due to the computation of 
possible collisions and cell crossing times scales as O(N) 
instead of 0(N 2 ), while the cost due to the event schedul- 
ing is 0(N log N) per unit of physical time instead of 
0(N 3 ) — a huge reduction. 

Contrary to what their scaling may suggest, one of- 
ten finds that the cost of the computation of collision 
times greatly dominates the scheduling cost for finite N. 
This is due to fact that the computations of many of the 
collision times requires numerical root searches, although 
some can, and should, be done analytically. Thus, to gain 
further computational improvements, one has to improve 
upon the efficiency of the numerical search for collision 
event times. A non-standard time-saving technique that 
we have developed for this purpose is to use virtual col- 
lision events. In this case, the grid search (see Sec. Illlf) 
for a possible collision time of atoms i and j is carried 
out only over a fixed small number of grid points, thus 
limiting the scope of the root search to a small search 
interval. If no collisions are detected in this search in- 
terval, a virtual collision event is scheduled in the binary 
event tree, much as if it were a possible future collision 
at the time of the last grid point that was investigated. 
If the point at which the grid search is curtailed is rather 
far in the future, it is likely this virtual event will not be 
executed because the atoms i and j probably will have 
collided with other atoms beforehand. Thus, computa- 
tional work has been saved by stopping the grid search 
after a few grid points. Every now and then however, 
atoms i and j will not have collided with other atoms 
at the time at which the grid search was stopped. In 
this case, the virtual collision event in the tree is exe- 
cuted, which entails continuing the root search from the 
point at which the search was previously truncated. The 
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continued search again may not find a root in a finite 
number of grid points and schedule another virtual col- 
lision, or it may now find a collision. In either case the 
new event is scheduled in the tree. This virtual collision 
technique avoids the unnecessary computation of a col- 
lision time that is so far in the future that it will not 
be executed in all likelihood anyway, while at the same 
time ensuring that if, despite the odds, that collision is to 
happen nonetheless, it is indeed found and correctly ex- 
ecuted. The trade-off of this technique is that the event 
tree is substantially larger, slowing down the event man- 
agement. Due to the high cost of numerical root searches 
however, the simulations presented in the accompanying 
paper showed that using virtual collision events yields an 
increase in efficiency between 25% and 110%, depending 
mainly on the system size. 



IV. COLLISION RULES 

At each moment of collision, the impulsive forces and 
torques lead to discontinuous jumps in the momenta and 
angular momenta of the colliding bodies. In the presence 
of constraints, there are two equivalent ways of deriv- 
ing the rules governing these changes. In the first ap- 
proach, the dynamics are treated by applying constraint 
conditions to Cartesian positions and momenta. This 
approach is entirely general and is suited for both con- 
strained rigid and non-rigid motion. In its generality, it is 
unnecessarily complicated for purely rigid systems and is 
not suitable for continuum bodies. The second approach, 
suitable for rigid bodies only, uses the fact that only six 
degrees of freedom, describing the center of mass motion 
and orientational dynamics are required to fully describe 
the dynamics of an arbitrary rigid body. The deriva- 
tion therefore consists of prescribing a collision process 
in terms of impulsive changes to the velocity of the center 
of mass and impulsive changes to the angular velocity. 



A. Constrained variable approach 

The general collision process in systems with discon- 
tinuous potentials can be seen as a limit of the collision 
process for continuous systems in which the interaction 
potential becomes infinitely steep. A useful starting point 
for deriving the collision rules is therefore to consider the 
effect of a force applied to the overall change in the mo- 
mentum of any atom k: 

p' k = Pk {t c + At) = Pk + / p k (t)dt 

Jt c -At 

pt e +At 

= Pk+ F k (r N {t))dt, (67) 

Jt c -At 

where Fk is the total force acting on atom k and pk = 
Pk{t c — At). Furthermore, here and below the pre and 



post-collision values of a quantity a are denoted by a and 
a', respectively. 

For discontinuous systems, the intermolecular forces 
are impulsive and occur only at an instantaneous colli- 
sion time t c . When atoms i and j collide, the interaction 
potential $ depends only on the scalar distance be- 
tween those atoms, so that the force on an arbitrary atom 
k is given by (without summation over i and j) 

d^(r N ) = d^(r N )dr lJ dr rj 
dr k dr i0 dr i: 

_ d$(r N ) 1 

drij r.ij 



dr k 



rij{Sjk - 8ik). 



Note that this is non-zero only for the atoms involved in 
the collision, as expected. Given that the force is impul- 
sive, it may be written as 



ggKrjy) 
drk 



— SS(t — t c )fij(Sjk - Sik) 7 



(68) 



where the scalar S is the magnitude of the impulse (to 
be determined) on atom a in the collision. 

In general, the constraint forces on the right-hand side 
of Eq. |J7Jl must also have an impulsive component when- 
ever intermolecular forces are instantaneous in order to 
maintain the rigid body constraints at all times. We ac- 
count for this by writing the Lagrange multipliers as 

K = v a + fj, a S(t - t c ). 

Because X a enters into the equations of motion for all 
atoms k involved in the constraint a a , there is an effect of 
this impulsive constraint force on all those atoms. Thus, 
one can write for the force on a atom k when atoms i and 
j collide: 



F k (r N ) = -v a 

+ S(t-t c ) 



da a (r N ) 
dr k 

Srij(Sjk - S ik ) - Ma 



da a (r N ) 
dr k 



(69) 



Substituting Eq. Ij69(l into l|67|) . one finds that the term 
proportional to v a vanishes in the limit that the time 
interval At approaches zero, so that the post-collision 
momenta p' k are related to the pre-collision momenta pk 
by 

Pk =Pk Ma dMi-jv) + Srij(5 jk - 5 ik ). (70) 

Note that at the instant of collision t — t c , the positions 
of all atoms r^r remain the same (only their momenta 
change) so that there is no ambiguity in the right-hand 
side of Eq. I|70l) as to whether to take the before or 
after the collision. It is straightforward to show that due 
to the symmetry of the interaction potential, the total 
linear momentum and angular momentum of the system 
are conserved by the collision rule Ea. l|70|l for arbitrary 
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values of the unknown scalar functions S and fi a . In ad- 
dition to these constants of the motion, the collision rule 
must also conserve total energy and preserve the con- 
straint conditions, a a {r^) = and a a {r^) = 0, before 
and after the collision. The first constraint condition is 
trivially satisfied at the collision time, since the positions 
are not altered at the moment of contact. The second 
constraint condition allows the scalar [i a to be related to 
the value of S using Eq. |j5Jl before and after the collision, 
since we must have 

^ = 2^— ■ TwT = ° 



For finite values of A<i>, the value of S is therefore 



m k dr k 
4^ m k dr k 



(71) 



-b± %/6 2 -4aA$ 
2a ' 



(77) 



where the physical solution corresponds to the positive 
(negative) root if b > (b < 0), provided b 2 > 4a A$. 
If this latter condition is not met, there is not enough 
kinetic energy to overcome the discontinuous barrier, 
and the system experiences a hard-core scattering, with 
A$ = 0, so that Eq. {ZZJ gives S = —b/a. Once the value 
of S has been computed, the discrete changes in momenta 
or velocities are easily computed using Eq. (|74|l . 



B. Rigid body approach 



Inserting Eq. (j7U|l into Eq. ij7T|) . one gets 

da 







V — (Srij(8jk -S ik )- ^/3^- 
^ rnfe V dr k 



Solving this linear equation for fi a gives 



= Sr 



1 dap 

m 3 dr 3 



1 da 



mi dr.; 



da a 
dr k 



.(72) 



(73) 



where the Z matrix was defined in Eq. (|12fl . Note that if 
atoms i and j are on different bodies, a given constraint 
ap involves either one or the other atom (or neither), so 
at least one of the two terms on the right-hand side of 
Eq. (SUl is then zero. Equation (fTTTf) can now be written 
as 



Pk 



p k + SAp k 



da 

Ap k = r^Sjk-Sik)-^-^, (74) 

dr k 

where /i* = fi a /S is a function of the phase-space coor- 
dinate as determined by Eq. (|73|l and is independent of 
S. 

Finally, the scalar S can be determined by employing 
energy conservation, 



The solution method outlined above can be applied to 
semi-flexible as well as rigid molecular systems, but is 
not very suitable for rigid, continuous bodies composed 
of an infinite number of point particles. For perfectly 
rigid molecules, a more convenient approach is therefore 
to analyze the effect of impulsive collisions on the center 
of mass and angular coordinates of the system, which are 
the minimum number of degrees of freedom required to 
specify the dynamics of rigid systems. The momentum 
of the center of mass P a and the angular momentum L a 
of rigid molecule a are affected by the impulsive collision 
via 



P' a = P a + AP a 

L' a = R a x P' a +I a -u}' a 

= L a + R a x AP a + I a ■ Auj a 

= L a + AL a , 



(78) 



where I a and u) a are the moment of inertia tensor and 
the angular velocity of body a in the laboratory frame, 
respectively. Note that they are related to their respec- 
tive quantities in the principal axis frame (body frame) 
via the matrix A. a (t) (now associated with the body a): 



T A T T A 

*-a — a a a 

U>a = A!w„. 



(79) 



2m k 



A$ 



Pk Pk 

2m k 



(75) 



where A$ = $' — $ denotes the discontinuous change in 
the potential energy at the collision time. Inserting the 
expression in (|74|) into l|75|l and using Eq. 10 , one gets a 
quadratic equation for the scalar S, 



aS 2 



bS + A$ = 
\ - Ap k ■ Ap k 



k 

E 



2m k 

Pk ■ Ap k 

m k 



(76) 



To derive specific forms for the impulsive changes 
AP a and Au> a , one may either calculate the impulsive 
force and torque acting on the center of mass and an- 
gular momentum, leading to AP a = — APf, = — Sr and 
AL a = —ALf, = r a x AP a , where r a and rj, are the 
points at which the forces are applied on body a and 6, 
respectively, while r = (rt,— r a )/|r(,— r a \, and S should be 
obtained from energy conservation. To understand this 
better and make a connection with the previous section, 
one may equivalently view the continuum rigid body as a 
limit of a non-continuum rigid body composed of n con- 
strained point particles, and use the expressions derived 
in the previous section for the changes in momenta of the 
constituents. 
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In the latter approach, it is convenient to switch the no- 
tation for the positions and momenta of the atoms from 
j-j and Pi to rf and pf, which indicate the position and 
momentum of particle i on body a, respectively. Using 
this notation and considering a collision between particle 
i on body a and particle j on body b, Eq. Ij74(l can be 
written as 

P t'-pt = SA P t - />a - 



-S 



rffS ik 



Mo 



dr k 



(80) 



1j is the unit vector along the direction 



where rfj = r"" jr 
of the vector rfj — — rf connecting atom i on body a 
with its colliding partner j on body b. Thus, noting that 

R a = Efc m t r k/ M a, where M a = Yl=i m k is the total 
mass of body a, and using Eq. JSDJ, one finds that 



P' a = P a + AP a 

n n 

AP a = Sj2&Pt = -Sv* a Yl 



(81) 



k=l 



k=l 



dcr a 
dr k 



Sfff 



-Sfff 



since 



E 

fc=i 



da a {r uv ) 



k=l 



drk 

Similarly, one finds that 
AL a 



cr^r^ 



(Sku - Skv) = o. 



(82) 



S^rfxApt 



k=l 



Ra X AP a 



S^r a k xAp k 

k=l 



where r\ 
evident that 



i- 



R„ 



x AP„ -Srf 



X TV 



(83) 



Comparing with Eq. I|78() . it is 



Au) n 



-si: 



(84) 



Note that I" 1 is a matrix inverse. Once again the im- 
pulsive changes are directly proportional to S, and the 
change of the angular velocity of body b in the laboratory 
frame due to the collision can be calculated analogously. 

To determine the scalar S, one again uses the conser- 
vation of total energy (E' — E) to see that 

IP'I 2 



\n\ 2 



2M a 2M b 



1 



1 



\P a \ 2 , IAI 2 i 



la U n 



1 

— U)h 

2 



2M a ' 2Mb ' 2" a - LaU}a+ 2" b - Ib " b - (85) 

Inserting Eqs. (|81|l and (|84|) into the energy equation 
above yields, after some manipulation, a quadratic equa- 
tion for S of the form of Eq. (|77|) , with 

AEt 



1 

2M a 

n .ab 



l 

2Mb 



AE a , 



where 



A E b , - nt ■ \l 



n 



with 



n 



n ij ~ r i x r< ^- 



For a spherically symmetric system, I = I = ill, and 



AE a ; = ri 



a,b\ 



n i / J l 



V. DYNAMICS IN THE CANONICAL AND 
MICRO-CANONICAL ENSEMBLES 

Any event-driven molecular dynamics simulation re- 
lies on the assumption that no collision is ever missed. 
However, collisions will be missed whenever the time dif- 
ference between two nearby events is on the order of 
(or smaller than) the time error of the scheduled events, 
which indicates that there is still a finite chance that a 
collision is missed even when event times are calculated 
in a simulation starting from an analytic expression, due 
to limits on machine precision. Although this subtle is- 
sue is not very important in a hard sphere system, in the 
present context it is of interest. Indeed, the extensive use 
of numerical root searches for the event time calculations 
combined with the need for computational efficiency de- 
mands a lower precision in the time values of collision 
events (typically a precision of 10 -10 instead of 10 -16 for 
analytical roots). In this section, it will be shown how 
to handle missed collisions in the context of the hybrid 
Monte Carlo scheme (HMC). 

In general, the HMC method 0, ^Tj combines 
the Monte Carlo method with molecular dynamics 
to construct a sequence of independent configurations 
{rjy , . . . , rjy }, distributed according to the canonical 
probability density 



P(fN) 



exp 



' fc R T 



(86) 



where Z is the configurational integral, fee is Boltzmann's 
constant, and T is the temperature. In the present con- 
text, this method can be implemented as follows: Ini- 
tially, a new set of momenta p' N is selected by choosing a 
random center of mass momentum P and angular veloc- 
ity u) for each molecule from the Gaussian distribution 



Pq oc exp 



1 



kuT \ 2M 



1 



-u-Iu 



The system is then evolved deterministically through 
phase-space for a fixed time r according to the 
equations of motion. This evolution defines a 
mapping of phase-space given by (rjv(0), Pjv(0)) *— > 



14 



{r N {To),pjf(ro)) = (r' N ,p' N ). The resulting phase space 
point (r' N ,p' N ) and trajectory segment are then accepted 
with probability 



PA(r' N ,p' N \r N (0),p N (0)) = min<{ l,exp 
where 



AE 



AE = E{r' N ,p' N ) - E(r N (Q),p N (0)), 



and 



A' 



E(r N ,p N ) = J2 ^rlP'l 2 + $ ( rAr ) 



i=l 



(87) 



(88) 



(89) 



This algorithm generates a Markov chain of configura- 
tions with an asymptotic distribution given by the sta- 
tionary canonical distribution defined in Eq. <|86|) pro- 
vided that the phase space trajectory is time reversible 
and area preserving |46||. Since free translational motion 
is time reversible, and the reversibility of the rotational 
equations of motion is evident from Eq. (|3()|) , the first re- 
quirement is satisfied. Furthermore, since the invariant 
phase space metric is uniform for fully rigid bodies (see 
Eq. (O in Sec. ITT"£|) . the area preserving condition is 
also satisfied. 

Ideally, a DMD simulation satisfies AE — so that ac- 
cording to Eq. l|S7|l every trajectory segment is accepted. 
In the less ideal, more realistic case in which collisions are 
occasionally missed, the HMC scheme provides a rigorous 
way of accepting or rejecting the segment. If a hard-core 
collision has been missed and the configuration at the end 
of a trajectory segment has molecules in unphysical re- 
gions of phase space where the potential energy is infinite, 
then AE = oo and the new configuration and trajectory 
segment are always rejected. On the other hand, if only 
a square-well interaction has been missed, AE at the end 
of the trajectory segment is finite and there is a non-zero 
probability of accepting the trajectory. 

An analogous strategy can be devised to carry out mi- 
crocanonical averages. In this case, the assignment of 
new initial velocities in the first step is still done ran- 
domly but in such a way that the total kinetic energy of 
the system remains constant. Such a procedure can be 
carried out by exchanging center of mass velocities be- 
tween randomly chosen pairs of molecules. The system 
is evolved dynamically through phase space for a fixed 
time To and the new phase space point (r' N ,p' N ) is ac- 
cepted according to 



Pa (r' N ,p' N \r N (0),p N (0)) 



if AE j£ 

1 if AE = 0. 



(90) 



where AE is given by Eq. (|S8|) . Clearly, the case AE ^ 
only occurs when a collision has been missed, and in such 
a case the trajectory segment is never accepted. 

It should be emphasized that in the HMC scheme, a 
new starting configuration for a segment of time evolu- 
tion is chosen only after every DMD time interval tq. An 



algorithm in which a new configuration is selected only 
after a collision is missed is likely to violate detailed bal- 
ance, and is therefore not a valid Monte-Carlo scheme. 
On the other hand, the length of the trajectory segments 
To in the HMC method outlined above can be chosen to 
be slightly larger than the relevant relaxation time of the 
system. Such a choice allows one to use the deterministic 
phase space trajectory to compute time-dependent cor- 
relation functions from the exact dynamics of the system 
without rejecting a significant fraction of the trajectory 
segments. 



VI. CONCLUSIONS 

In this paper we have shown how to carry out dis- 
continuous molecular dynamics simulations for arbitrary 
semi- flexible and rigid molecules. For semi-flexible bod- 
ies, the dynamics and collision rules have been derived 
from the principles of constrained Lagrangian mechan- 
ics. The implementation of an efficient DMD method 
for semi-flexible systems is hindered by the fact that in 
almost all cases the equations of motion must be prop- 
agated numerically in an event searching algorithm so 
that the constraints are enforced at all times. Nonethe- 
less, such a scheme can be realized using the SHAKE pgj 
or RATTLE 48j algorithms in combination with the root 
searching methods outlined here. 

The dynamics of a system of completely rigid molecules 
interacting through discontinuous potentials is more 
straightforward. For such a system, the Euler equations 
for rigid body dynamics can be used to calculate the free 
evolution of a general rigid object. This analytical solu- 
tion enables the design of efficient numerical algorithms 
for the search for collision events. In addition, the col- 
lision rules for calculating the discontinuous changes in 
the components of the center of mass velocity and angular 
momenta have been obtained for arbitrary bodies inter- 
acting through a point based on conservation principles. 
Furthermore, the sampling of the canonical and micro- 
canonical ensembles, as well as the handling of missed 
collisions, has also been discussed in the context of a hy- 
brid Monte Carlo scheme. 

From an operational standpoint, the difference be- 
tween the method of DMD and molecular dynamics us- 
ing continuous potentials in rigid systems lies in the fact 
that the DMD approach does not require the calculation 
of forces and sequential updating of phase space coor- 
dinates at discrete (and short) time intervals since the 
response of the system to an impulse can be computed 
analytically. Instead, the computational effort focuses 
on finding the precise time at which such impulses exert 
their influence. The basic building block outlined here 
for the numerical computation of collision times is a grid 
search, for which the positions of colliding atoms on a 
given pair of molecules need to be computed at equally 
spaced points in time. As outlined in Sec. II I II this can 
be done efficiently starting with a completely explicit an- 
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alytical form of the motion of a torque-free rigid body, 
without which the equations of motions would have to 
be integrated numerically. An efficient implementation 
of the DMD technique to find the time collision events 
should make use of a) a large grid step combined with 
a threshold scenario to catch pathological cases, b) so- 
phisticated but standard techniques such as binary event 
trees, cell divisions, and local clocks, and c) a new tech- 
nique of finding collision times numerically that involves 
truncating the grid search and scheduling virtual collision 
events. 

On a fundamental level, it is natural to wonder whether 
the 'stepped' form of a discontinuous potential could pos- 
sibly model any realistic interaction. Such concerns are 
essentially academic, since it is always possible to approx- 
imate a given interaction potential with as many (small) 
steps as one would like in order to approximate a given 
potential to any desired level of accuracy[|j. Of course, 
the drawback to mimicking a smooth potential with a 
discontinuous one with many steps is that the number of 
'collision' events that occur in the system per unit time 



scales with the number of steps in the potential. Hence, 
one would expect that the efficiency of the simulation 
scales roughly inversely with the number of steps in the 
interaction potential. Nonetheless, the issue is a practical 
one: How small can the number of steps in the interaction 
potential be such that one still gets a good description 
of the ph ysic s under investigation? In the accompany- 
ing paper |49j, we will see for benzene and methane that 
it takes surprisingly few steps (e.g. a hard core plus a 
square-well interaction) to get results which are very close 
to those of continuous molecular dynamics. Additionally, 
we compare the efficiency of such simulations to simula- 
tions based on standard molecular dynamics methods. 
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